Effects of pressure on diffusion and vacancy formation in MgO from non-empirical 

free-energy integrations. 
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The free energies of vacancy pair formation and migration in MgO were computed via molecular 
dynamics using free-energy integrations and a non-empirical ionic model with no adjustable param- 
eters. The intrinsic diffusion constant for MgO was obtained at pressures from to 140 GPa and 
temperatures from 1000 to 5000 K. Excellent agreement was found with the zero pressure diffusion 
data within experimental error. The homologous temperature model which relates diffusion to the 
melting curve describes well our high pressure results within our theoretical framework. 
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Diffusion and vacancy formation are critical to kinetic 
processes in materials, yet little is known about diffusion 
at ultra-high pressures due to experimental difficulties. 
Rheology of oxide minerals at high pressures is also cru- 
cial in geophysics and is dependent on diffusive behav- 
ior which is only available experimentally at relatively 
low pressures [0. In ionic systems such as MgO, the 
dominant intrinsic defect is the pair vacancy |2|-||] with 
Mg and O sites vacant. Previous work on vacancies in 
MgO used pseudopotential computations or lattice dy- 
namics or the Mott-Littleton approach with a variety of 
semi- empirical potentials . The accuracy of quasi- 
harmonic lattice dynamics calculations degrades above 
the Debye temperature and the Mott-Littleton procedure 
and pseudopotential computations were restricted to K. 

We used molecular dynamics (MD) with non-empirical 
potentials to determine the self-diffusion coefficient D 
where [0 



D 



6 



l v exp( 



AG/ 

w 



AG,, 



k b T 



(1) 



Zf is the number of equivalent ways of forming a vacancy 
type, Z m is the number of equivalent diffusion paths, I 
is the jump distance, v is the attempt frequency, AG/ 
and AG m are the energies of formation and migration, 
respectively and W is the solubility factor for polyatomic 
materials. If the sites are uncorrelated (Schottky defect), 
then, for rocksalt structured (Bl) crystals such as MgO, 
W = 2, Zf = 1, Z rn = 12. Highly correlated defects 
(bound pair) require W =1 and Zt = 6. Symmetry and 
energy considerations determine the value of Z m . In ei- 
ther case, l 2 = a 2 /2 where a is the cubic cell parameter. 

We used the variational induced breathing (VIB) 
model which reliably gives the thermal properties and 
equation of state of MgO || to compute the energetics 
and interatomic forces. The VIB model is a Gordon- 
Kim type model H in which the total charge density is 



modeled by overlapping ionic charge densities which are 
computed using the local density approximation (LDA) 
|i"of . The total energy is a sum of three terms: (a) 
the long-range electrostatic energy computed using the 
Ewald method, (b) the self-energy of each atom and (c) 
the short range interaction energy, the sum of the ki- 
netic, short-range electrostatic and exchange-correlation 
energies from the LDA. There are three approximations 
beyond the LDA (1) The charge density is modeled 
rather than computed self-consistently. Comparisons 
with accurate linearized augmented plane wave (LAPW) 
computations show this is a good approximation for MgO 
]i"2| ; (2) the pair approximation is used for the short- 
range interactions (c) which is a good approximation as 
long as closed shell ions are used Q; (3) the Thomas- 
Fermi kinetic energy is used for the short-range overlap 
kinetic energy. The self-energy (b) includes the correct 
LDA Kohn-Sham kinetic energy. O 2- is not stable in the 
free state and is stabilized by introducing a sphere of 2+ 
charge (Watson sphere) around it in the LDA atomic cal- 
culations. Interactions are obtained for overlapping ion 
pairs at different distances with different Watson sphere 
radii on the O's. For efficiency, the interactions were fit 
with a 21 parameter analytical expression as functions of 
r, the interatomic distance, and E/j = Zi/Ri where Ui,Zi, 
and Ri are the Watson sphere potential, charge (2+) and 
radius for atom i, respectively. During the simulations, 
the total energy was variationally optimized with respect 
to all of the Watson sphere radii at each time step. 

The attempt frequency v was determined by Fourier 
transforming the trajectories of the diffusing ion pro- 
jected onto the shortest path to the vacancy. We con- 
sidered two models: (1) the lowest frequency peak in the 
spectrum, assuming that the diffusive motion is mostly 
from the lowest energy mode and (2) the average fre- 
quency computed from the Fourier transform. In the 
first case we found that v = 5.2 THz and is independent 
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of P and T over the range studied. The second case gave 
an attempt frequency that is within a factor 2 of the low 
frequency value . Given the uncertainties in the calcula- 
tions and experimental determinations, the difference in 
the final results between these two approaches is small 
and we adopted case (1) below. 

Free energies were computed with the finite time vari- 
ational or "adiabatic switching" thermodynamic integra- 
tion method [E3 . The free energy difference between the 
initial and final state is 
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where A is a progress variable which ranges from to 1 as 
the system "switches" from its initial to final state, H is 
the system Hamiltonian, and (), represents an ensemble 
average. In order to obtain AG/, we first calculated the 
free energy difference between an ideal crystal at volume, 
Vi, giving the desired average P at T and an Einstein 
crystal at the same V and T. This was repeated for 
a defective crystal with a bound vacancy pair in each 
periodic cell at Vd corresponding to P. Then for an N 
atom periodic cell 
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where Fp~ L (Vd) is the Helmholtz free energy of a de- 
fective crystal with L vacant sites and Fj (Vj) is the 
Helmholtz free energy of the ideal crystal. The Hamilto- 
nian took the form 



H(X)=H Y1B x (l-\) + H ein \. 



(4) 



where H e i n is the Hamiltonian for an Einstein crystal Jll 
which can be written as 
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where K is the kinetic energy, U a is the static contribu- 
tion to the potential, mi, x\, x^o, and uj ein i are the mass, 
position, static lattice position, and Einstein frequency of 
the zth particle, respectively. The form of A as a function 
of the scaled time, r, is 

A(r) = r 5 (70r 4 - 315r 3 + 540r 2 - 420r + 126) (6) 

where r = t/t s , t is the elapsed time, and t s is the total 
switching time |15| . 

Migration free energies were calculated using the adi- 
abatic switching procedure at constant P and T. We 
computed the energy it takes to push the atom out of 
one lattice site and into another vacant lattice site Jl6| . 
The force on the migrating atom due to Hyi b in the mi- 
gration direction was set to zero and the negative of this 
force was evenly distributed among the rest of the atoms 



so that the force on the center of mass was zero. The 
position of the migrating atom was then incremented in 
the migration direction. Forces on the the atom in the 
plane perpendicular to the migration direction were not 
artificially constrained so that the path the migrating 
atom took did not lie on a direct line between the initial 
position of the atom and the vacancy site. AG m was de- 
termined by summing the difference in -ffyiB before and 
after incrementing the position of the migrating atom 
with the other atomic positions held fixed. We found 
that the barriers to migration for ions that do not have 
a vacancy as a nearest neighbor were lower than for ones 
who do. Thus the value of Z m is 8 in MgO. 

MD was performed using a timestep of 1 fs with a 5th 
order Gear predictor-corrector scheme JT^ in an isobaric- 
isothermal ensemble generated using the extended sys- 
tem method Ul for 10 ps equilibration times followed by 
a 10(15) ps switching time for the formation(migration) 
energy. Convergence with respect to our nominal 216 
atom system size was verified for systems with up to 
1000 atoms to be within 1%. Doubling the integra- 
tion time resulted in free energy variations of 1% while 
halving it increased the calculated free energy change 
by 10%. The computationally efficient first principles 
method used here lends itself to demanding convergence 
tests, especially with respect to system size, that would 
prove too time-consuming with self-consistent methods 

Values of AG/ for bound pairs and AG m are given 
in Table Q. At GPa, these energies are within 5% of 
those derived from previous theoretical and experimen- 
tal results |^f!||- To determine the dominant vacancy 
mechanism, we calculated the binding energy AG;, of a 
bound pair from the difference in AG/ between bound 
vacancy pairs and those with the largest possible dis- 
tance between vacant Mg and O sites in a 1000 atom 
supercell corrected for image forces [^()| and found AG& 
= 2.5xlO~ 19 J at GPa and 2000 K and 7.0xlO~ 19 J at 
140 GPa and 3000 K. No significant changes in migration 
energy relative to the bound pair simulation were found. 
Assuming that the variation in binding energy is linear 
in pressure and independent of temperature and taking 
into account the configurational entropy, we calculated 
the vacancy concentration and Gibbs free energy change 
for crystals containing bound and disassociated pairs rel- 
ative to the perfect crystal M. This analysis shows that 
Schottky defects dominate at 1000 K or above due to 
entropy contributions (lower temperatures will favor the 
bound state). Ionic conductivity measurements indicate 
that Mg diffusion (DMg) is controlled by impurities |l9| ] 
whereas O diffusion (Do) is intrinsic in nature and di- 
rectly comparable to our analysis. We find that our pre- 
dicted Do agrees with experiment within mutual error 
(Fig- 1). 

Fitting our diffusion constants with the relation 
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TABLE I. Free Energies of Formation and Migration and 
Cell Parameter 



20 40 60 80 100 120 140 
Pressure (GPa) 

FIG. 1. Predicted pressure and temperature dependence 
of the self-diffusion coefficients in MgO. Curves represent the 
best fit to the coefficients using the activation energy-volume 
relation given by Equation ^. Vertical symbol size of experi- 
mental datum taken from Ref. W represents uncertainty. 



hxD = \n{a 2 v) + S* +PS* ' 
-(E* + PV Q * + P 2 V *')/k b T 
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gives the zero pressure activation entropy S* =3- 
(4)fc(,, its pressure derivative S*' =0.03-(0.02)fc;„ activa- 
tion energy E* =9.0-(9.4) xl0~ 19 J, activation volume 
V* = 16.0-(16.7) A 3 , and its pressure derivative V*' =- 
0.031- (-0.038) A 3 /GPa for Mg-(O). The activation vol- 
ume varies as a function of pressure consistant with pre- 
vious discussions ]2C|-[2^]. The activation entropy is of 
the same order as previous estimates of the formation 
entropy at GPa but varies much less drastically with 
pressure than a previous estimate Q . 

Finally, we considered the homologous temperature re- 
lation 



D = D exp(gT m /T) 



(8) 



commonly used to model the dependence of diffusion on 
P and T J2j|. Because of the similarity of behavior in 
diffusion of Mg and O, we used the effective diffusion co- 
efficient, D eff = 2D Mg D /(D Mg + D ) for D @. We 
tested this model using the theoretical melting curve of 
Cohen and Weitz j25| obtained with the same VIB po- 
tential as used here, and the extrapolated experimental 
melting curve of Zerr and Boehler |p6| which has a lower 
dP/dT. Good global fits were found using the theoret- 
ical melting curve, but the experimental melting curve 
is not consistent with the present diffusion results. As 
discussed in Ref. pq| the experimental results may be in- 
fluenced by Ar solubility in MgO melt at high pressures. 



p 


T 


AG/ 


a 


AG 




(GPa) 


(K) 




(A) 


(10~ 19 J) 










Mg 


O 


o 


1000 


8 1 9 + 39 


4.2495 


3.36 ± 0.10 


3.75 ± 0.11 





2000 


7.88 ± 0.23 


4.3073 


2.72 ± 0.12 


3.16 ± 0.12 


20 


2000 


12.04 ± 0.40 


4.1312 


3.75 ± 0.10 


4.09 ± 0.10 


20 


3000 


11.57 ± 0.56 


4.1687 


3.60 ± 0.10 


4.06 ± 0.12 


80 


2000 


21.20 ± 0.60 


3.8610 


5.16 ± 0.10 


5.84 ± 0.17 


80 


3000 


21.01 ± 0.90 


3.8803 


4.99 ± 0.11 


5.86 ± 0.07 


80 


4000 


20.75 ± 0.71 


3.9031 


4.55 ± 0.10 


5.61 ± 0.16 


140 


2000 


26.96 ± 0.82 


3.7071 


6.39 ± 0.10 


6.91 ± 0.10 


140 


3000 


25.29 ± 1.11 


3.7210 


5.99 ± 0.15 


6.21 ± 0.21 


140 


5000 


23.44 ± 1.76 


3.7509 


4.80 ± 0.30 
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Theoretical estimates of the melting curve are generally 
consistent with each other and with expected thermody- 
namic parameters. 

We also tested the use of zero pressure diffusion re- 
sults only in Eq. || and found that extrapolations to high 
pressure using the melting curve were reasonably reliable 
although some accuracy was lost compared to the results 
from direct high pressure simulations. This gives justifi- 
cation for use of melting curves in estimating high P and 
T diffusion in oxides. In addition, we found that g was 
less than 14 while the average value for alkali halides is 
24 [^7| indicating that conclusions based on a systematic 
value for this parameter may be invalid. 

In summary, we found (1) excellent agreement with 
experimental results, (2) that defects are formed from 
Schottky pairs as opposed to neutral divacancies, and 
(3) the homologous temperature relation holds within our 
theoretical framework. These results will help constrain 
rheological properties of the deep Earth and provide con- 
straints for pressure effects on kinetics in oxides. 
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